Two independent samples T test
·
Based on: Handbook of Parametric and Nonparametric Statistical Procedures, 5th edition, Chapter 11
using Random
Random.seed!(783376894512)
using Distributions
using Intervals
using Plots
using LinearAlgebra
n_sim = 1_000_000 # number of simulations to verify theoretical distributions
True model
n1 = 4
n2 = 5
σₜ = 2.0
μₜ1 = 3.0
μₜ2 = 4.0
Mₜ1 = Product([Normal(μₜ1, σₜ) for _ in 1:n1]) # assumed known to be Normal and i.i.d.
Mₜ2 = Product([Normal(μₜ2, σₜ) for _ in 1:n2]) # assumed known to be Normal and i.i.d.
generate_data(M1, M2) = (rand(M1), rand(M2))
Hypothesis testing
function test_statistic(y1, y2, Δμ₀)
n1 = length(y1)
n2 = length(y2)
μ̂1 = sum(y1) / n1
μ̂2 = sum(y2) / n2
σ̂²1 = sum((y1 .- μ̂1) .^ 2) / (n1 - 1)
σ̂²2 = sum((y2 .- μ̂2) .^ 2) / (n2 - 1)
σ̂²p = ((n1 - 1) * σ̂²1 + (n2 - 1) * σ̂²2) / (n1 + n2 - 2)
return (μ̂1 - μ̂2 - Δμ₀) / (sqrt(σ̂²p) * sqrt(1 / n1 + 1 / n2))
end
Null Hypothesis is true
Δμ₀ = -1.0
Δμ₀ == μₜ1 - μₜ2
all(mean(Mₜ1) .== μₜ1)
all(mean(Mₜ2) .== μₜ2)
all(sqrt.(var((Mₜ1))) .== σₜ)
all(sqrt.(var((Mₜ2))) .== σₜ)
simulated_statistics = [test_statistic(generate_data(Mₜ1, Mₜ2)..., Δμ₀) for _ in 1:n_sim]
distribution_under_H₀ = t -> pdf(TDist(n1 + n2 - 2), t)
Code
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n when Null Hypothesis is correct",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Null Hypothesis is true, but assumptions don’t hold
Observations are non-normal
M1_non_normal = Product([
MixtureModel(
[Normal(μₜ1 + sqrt(3.99), 0.1), Normal(μₜ1 - sqrt(3.99), 0.1)], [0.5, 0.5]
) for _ in 1:n1
])
M2_non_normal = Product([
MixtureModel(
[Normal(μₜ2 + sqrt(3.99), 0.1), Normal(μₜ2 - sqrt(3.99), 0.1)], [0.5, 0.5]
) for _ in 1:n2
])
all(mean(M1_non_normal) .== μₜ1)
all(mean(M2_non_normal) .== μₜ2)
all(sqrt.(var(M1_non_normal)) .== σₜ)
all(sqrt.(var(M2_non_normal)) .== σₜ)
simulated_statistics = [
test_statistic(generate_data(M1_non_normal, M2_non_normal)..., Δμ₀) for _ in 1:n_sim
]
Code
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n under H₀, but non-normal observations",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Observations are non-independent
M1_non_independent = MvNormal(
[μₜ1, μₜ1, μₜ1, μₜ1, μₜ1],
[
σₜ^2 0.5*σₜ^2 0 0 0
0.5*σₜ^2 σₜ^2 0 0 0
0 0 σₜ^2 0 0
0 0 0 σₜ^2 0
0 0 0 0 σₜ^2
],
)
Σ2 = collect(Diagonal(fill(σₜ^2, n2)))
Σ2[2, 1] = Σ2[1, 2] = 0.5 * σₜ^2
M2_non_independent = MvNormal(fill(μₜ2, n2), Σ2)
all(mean(M1_non_independent) .== μₜ1)
all(mean(M2_non_independent) .== μₜ2)
all(sqrt.(var(M1_non_independent)) .== σₜ)
all(sqrt.(var(M2_non_independent)) .== σₜ)
simulated_statistics = [
test_statistic(generate_data(M2_non_independent, M2_non_independent)..., Δμ₀) for
_ in 1:n_sim
]
Code
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n under H₀, but non-independent observations",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Observations are non-identical
M1_non_identical = Product([Normal(μₜ1, i * σₜ) for i in 1:n1])
M2_non_identical = Product([Normal(μₜ2, i * σₜ) for i in 1:n1]) # note the n1 here to differentiate from the next counterexample
all(mean(M1_non_identical) .== μₜ1)
all(mean(M2_non_identical) .== μₜ2)
any(sqrt.(var(M1_non_identical)) .== σₜ)
any(sqrt.(var(M2_non_identical)) .== σₜ)
simulated_statistics = [
test_statistic(generate_data(M1_non_identical, M2_non_identical)..., Δμ₀) for
_ in 1:n_sim
]
Code
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n under H₀, but non-identical observations",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Two samples have a different variance
M2_different_variance = Product([Normal(μₜ1, 2 * σₜ) for _ in 1:n2])
all(mean(M2_different_variance) .== μₜ1)
all(sqrt.(var(M2_different_variance)) .!= σₜ)
simulated_statistics = [
test_statistic(generate_data(Mₜ1, M2_different_variance)..., Δμ₀) for _ in 1:n_sim
]
Code
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n under H₀, but different sample variances",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Null Hypothesis is false
Δμ₀ = 2.0
Δμ₀ != μₜ1 - μₜ2
simulated_statistics = [test_statistic(generate_data(Mₜ1, Mₜ2)..., Δμ₀) for _ in 1:n_sim]
distribution_under_H₁ =
tnc -> pdf(NoncentralT(n1 + n2 - 2, (μₜ1 - μₜ2 - Δμ₀) / √(σₜ^2 / n1 + σₜ^2 / n2)), tnc)
Code
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n when Null Hypothesis is incorrect",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Null Hypothesis is false, but assumptions don’t hold
Same counter examples as above.
Code
simulated_statistics = [
test_statistic(generate_data(M1_non_normal, M2_non_normal)..., Δμ₀) for _ in 1:n_sim
]
1000000-element Vector{Float64}:
-2.7319646451633637
-1.1635565113397672
-2.5973138575913843
-3.1753393673498325
-0.6945610789746434
-0.5825206489634959
-0.9541786713792212
-0.6833793462505765
-1.6870049418296202
-1.1047206335419304
⋮
-1.0827951572184749
-2.4820903467430555
-1.342302654041741
-3.0872701991748603
-1.0679875565707986
-3.2196581992897255
-0.7061209569844831
-1.6155670827847624
-0.580786774694805
Code
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n under H₁, but non-normal observations",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Code
simulated_statistics = [
test_statistic(generate_data(M2_non_independent, M2_non_independent)..., Δμ₀) for
_ in 1:n_sim
]
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n under H₁, but non-independent observations",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Code
simulated_statistics = [
test_statistic(generate_data(M1_non_identical, M2_non_identical)..., Δμ₀) for
_ in 1:n_sim
]
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n under H₁, but non-identical observations",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Code
simulated_statistics = [
test_statistic(generate_data(Mₜ1, M2_different_variance)..., Δμ₀) for _ in 1:n_sim
]
histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
title="Distribution of two independent samples T test statistic, \n under H₁, but different sample variances",
xlimits=(tmin, tmax),
xlabel="t",
ylabel="pdf",
grid=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Confidence Interval
α = 0.05 # proportion of CIs which should fail to cover true values
function confidence_interval(y1, y2, α)
n1 = length(y1)
n2 = length(y2)
μ̂1 = sum(y1) / n1
μ̂2 = sum(y2) / n2
σ̂²1 = sum((y1 .- μ̂1) .^ 2) / (n1 - 1)
σ̂²2 = sum((y2 .- μ̂2) .^ 2) / (n2 - 1)
σ̂²p = ((n1 - 1) * σ̂²1 + (n2 - 1) * σ̂²2) / (n1 + n2 - 2)
L = μ̂1 - μ̂2 + quantile(TDist(n1 + n2 - 2), α / 2) * (√(σ̂²p * (1 / n1 + 1 / n2)))
U = μ̂1 - μ̂2 + quantile(TDist(n1 + n2 - 2), 1 - α / 2) * (√(σ̂²p * (1 / n1 + 1 / n2)))
return Interval(L, U)
end
simulated_CIs = [confidence_interval(generate_data(Mₜ1, Mₜ2)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim
Code
bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test",
ylimits=(0.0, 1.0),
yticks=[0.05, 0.95],
ylabel="proportion",
grid=true,
legend=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
assumptions don’t hold
Same counter examples as above.
Code
simulated_CIs = [confidence_interval(generate_data(M1_non_normal, M2_non_normal)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim
bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test \n but non-normal observations",
ylimits=(0.0, 1.0),
yticks=[0.05, 0.95],
ylabel="proportion",
grid=true,
legend=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
#top_margin = 25Plots.px,
)
Code
simulated_CIs = [confidence_interval(generate_data(M1_non_independent, M2_non_independent)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim
bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test \n but non-independent observations",
ylimits=(0.0, 1.0),
yticks=[0.05, 0.95],
ylabel="proportion",
grid=true,
legend=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Code
simulated_CIs = [confidence_interval(generate_data(M1_non_identical, M2_non_identical)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim
bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test \n but non-identical observations",
ylimits=(0.0, 1.0),
yticks=[0.05, 0.95],
ylabel="proportion",
grid=true,
legend=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)
Code
simulated_CIs = [confidence_interval(generate_data(Mₜ1, M2_different_variance)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim
bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test \n but different sample variances",
ylimits=(0.0, 1.0),
yticks=[0.05, 0.95],
ylabel="proportion",
grid=true,
legend=false,
legendfontsize=12,
tickfontsize=12,
guidefontsize=16,
)